/* SPDX-License-Identifier: SunMicrosystems */
/* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. */

/**
 *
 * This family of functions implements the natural logarithm of :math:`1 + x`.
 *
 * Synopsis
 * ========
 *
 * .. code-block:: c
 *
 *     #include <math.h>
 *     float log1pf(float x);
 *     double log1p(double x);
 *     long double logl(long double x);
 *
 * Description
 * ===========
 *
 * ``log1p`` computes the natural logarithm of :math:`1 + x`.
 *
 * Mathematical Function
 * =====================
 *
 * .. math::
 *
 *    log1p(x) \approx ln(1 + x)
 *
 * Returns
 * =======
 *
 * ``log1p`` returns the natural logarithm of :math:`1 + x`.
 *
 * Exceptions
 * ==========
 *
 * Raise ``invalid operation`` exception when the input value is less than :math:`-1`.
 *
 * Raise ``divide by zero`` exception when the input value is :math:`-1`.
 *
 * .. May raise ``underflow`` exception.
 *
 * Output map
 * ==========
 *
 * +---------------------+---------------+---------------+---------------+-------------------+---------------+---------------+-------------------+---------------+---------------+
 * | **x**               | :math:`-Inf`  | :math:`<-1`   | :math:`-1`    | :math:`]-1,0[`    | :math:`-0`    | :math:`+0`    | :math:`>0`        | :math:`+Inf`  | :math:`NaN`   |
 * +=====================+===============+===============+===============+===================+===============+===============+===================+===============+===============+
 * | **log1p(x)**        | :math:`qNaN`  | :math:`qNaN`  | :math:`-Inf`  | :math:`ln(1 + x)` | :math:`x`                     | :math:`ln(1 + x)` | :math:`+Inf`  | :math:`qNaN`  |
 * +---------------------+---------------+---------------+---------------+-------------------+---------------+---------------+-------------------+---------------+---------------+
 *
 *///

/* double log1p(double x)
 *
 * Method :
 *   1. Argument Reduction: find k and f such that
 *            1+x = 2^k * (1+f),
 *       where  sqrt(2)/2 < 1+f < sqrt(2) .
 *
 *      Note. If k=0, then f=x is exact. However, if k!=0, then f
 *    may not be representable exactly. In that case, a correction
 *    term is need. Let u=1+x rounded. Let c = (1+x)-u, then
 *    log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u),
 *    and add back the correction term c/u.
 *    (Note: when x > 2**53, one can simply return log(x))
 *
 *   2. Approximation of log1p(f).
 *    Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
 *         = 2s + 2/3 s**3 + 2/5 s**5 + .....,
 *              = 2s + s*R
 *      We use a special Remez algorithm on [0,0.1716] to generate
 *     a polynomial of degree 14 to approximate R The maximum error
 *    of this polynomial approximation is bounded by 2**-58.45. In
 *    other words,
 *                2      4      6      8      10      12      14
 *        R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s  +Lp6*s  +Lp7*s
 *      (the values of Lp1 to Lp7 are listed in the program)
 *    and
 *        |      2          14          |     -58.45
 *        | Lp1*s +...+Lp7*s    -  R(z) | <= 2
 *        |                             |
 *    Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
 *    In order to guarantee error in log below 1ulp, we compute log
 *    by
 *        log1p(f) = f - (hfsq - s*(hfsq+R)).
 *
 *    3. Finally, log1p(x) = k*ln2 + log1p(f).
 *                  = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
 *       Here ln2 is split into two floating point number:
 *            ln2_hi + ln2_lo,
 *       where n*ln2_hi is always exact for |n| < 2000.
 *
 * Special cases:
 *    log1p(x) is NaN with signal if x < -1 (including -INF) ;
 *    log1p(+INF) is +INF; log1p(-1) is -INF with signal;
 *    log1p(NaN) is that NaN with no signal.
 *
 * Accuracy:
 *    according to an error analysis, the error is always less than
 *    1 ulp (unit in the last place).
 *
 * Constants:
 * The hexadecimal values are the intended ones for the following
 * constants. The decimal values may be used, provided that the
 * compiler will convert from decimal to binary accurately enough
 * to produce the hexadecimal values shown.
 *
 * Note: Assuming log() return accurate answer, the following
 *      algorithm can be used to compute log1p(x) to within a few ULP:
 *
 *        u = 1+x;
 *        if(u==1.0) return x ; else
 *               return log(u)*(x/(u-1.0));
 *
 *     See HP-15C Advanced Functions Handbook, p.193.
 */

#include <math.h>
#include "../common/tools.h"
#include "internal/log1pmfd.h"

#ifndef __LIBMCS_DOUBLE_IS_32BITS

static const double
ln2_hi  =  6.93147180369123816490e-01,  /* 3fe62e42 fee00000 */
ln2_lo  =  1.90821492927058770002e-10;  /* 3dea39ef 35793c76 */

static const double zero = 0.0;

double log1p(double x)
{
#ifdef __LIBMCS_FPU_DAZ
    x *= __volatile_one;
#endif /* defined(__LIBMCS_FPU_DAZ) */

    double hfsq, f, c, R, u;
    int32_t k, hx, hu, ax;

    GET_HIGH_WORD(hx, x);
    ax = hx & 0x7fffffff;

    k = 1;

    if (hx >= 0x7ff00000) {                /* x = NaN/+-Inf */
        return x + x;
    }

    if (hx < 0x3FDA827A) {                 /* x < 0.41422  */
        if (ax >= 0x3ff00000) {            /* x <= -1.0 */
            if (isnan(x)) {
                return x + x;
            } else if (x == -1.0) {
                return __raise_div_by_zero(-1.0); /* log1p(-1)=-inf */
            } else {
                return __raise_invalid(); /* log1p(x<-1)=NaN */
            }
        }

        if (ax < 0x3e200000) {             /* |x| < 2**-29 */
            if (ax < 0x3c900000) {         /* |x| < 2**-54 */
                return __raise_inexact(x);
            } else {
                return __raise_inexact(x - x * x * 0.5);
            }
        }

        if (hx > 0 || hx <= ((int32_t)0xbfd2bec3U)) {
            k = 0;
            f = x;
            hu = 1;
        }                                  /* -0.2929<x<0.41422 */
    }

    if (k != 0) {
        if (hx < 0x43400000) {
            u  = 1.0 + x;
            GET_HIGH_WORD(hu, u);
            k  = (hu >> 20) - 1023;
            c  = (k > 0) ? 1.0 - (u - x) : x - (u - 1.0); /* correction term */
            c /= u;
        } else {
            u  = x;
            GET_HIGH_WORD(hu, u);
            k  = (hu >> 20) - 1023;
            c  = 0;
        }

        hu &= 0x000fffff;

        if (hu < 0x6a09e) {
            SET_HIGH_WORD(u, hu | 0x3ff00000); /* normalize u */
        } else {
            k += 1;
            SET_HIGH_WORD(u, hu | 0x3fe00000); /* normalize u/2 */
            hu = (0x00100000 - hu) >> 2;
        }

        f = u - 1.0;
    }

    hfsq = 0.5 * f * f;

    if (hu == 0) { /* |f| < 2**-20 */
        if (f == zero) {
            c += k * ln2_lo;
            return k * ln2_hi + c;
        }

        R = hfsq * (1.0 - 0.66666666666666666 * f);

        return k * ln2_hi - ((R - (k * ln2_lo + c)) - f);
    }

    if (k == 0) {
        return f - (hfsq - __log1pmf(f));
    } else {
        return k * ln2_hi - ((hfsq - (__log1pmf(f) + (k * ln2_lo + c))) - f);
    }
}

#ifdef __LIBMCS_LONG_DOUBLE_IS_64BITS

long double log1pl(long double x)
{
    return (long double) log1p((double) x);
}

#endif /* #ifdef __LIBMCS_LONG_DOUBLE_IS_64BITS */
#endif /* #ifndef __LIBMCS_DOUBLE_IS_32BITS */
